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Abstract 

NEIGHBORING  EXTREMAL  GUIDANCE  FOR 
SYSTEMS  WITH  PIECEWISE  LINEAR 
CONTROL  USING  TIME  AS  THE 
REFERENCE  VARIABLE 

by 

WILLIAM  ABRAM  LIBBY,  B.S.A.E. 

SUPERVISING  PROFESSOR:  Dr.  David  G.  Hull 

A  guidance  law  for  the  control  of  a  system  in  the  neighborhood  of  a 
nominal  suboptimal  trajectory  is  developed.  The  guidance  law  is  demonstrated 
using  a  lunar  launch  problem  with  constraints  at  orbit  entry.  A  set  of 
precomputed  gains  is  used  by  the  guidance  law  to  operate  on  an  extremal  path  in 
the  neighborhood  of  the  suboptimal  trajectory.  The  guidance  law  and  gains  are 
designed  to  minimize  the  change  in  the  desired  performance  index  while  still 
satisfying  the  final  path  constraints. 

In  the  lunar  launch  problem,  the  nominal  suboptimal  trajectory  minimizes 
the  final  time  using  piecewise  linear  control.  This  trajectory  is  obtained  to 
provide  a  nominal  control  history.  The  guidance  law  is  found  by  minimizing  the 


second  variation  of  the  suboptimal  trajectory  performance  index  subject  to  the 
final  constraints  being  satisfied.  For  the  lunar  launch  problem,  the  guidance  law 
leads  to  a  set  of  gains  that  relates  deviations  from  the  suboptimal  trajectory  to 
required  changes  in  the  nominal  control  history.  The  deviations  from  the 
suboptimal  trajectory,  used  together  with  the  precomputed  gains,  determines  the 
change  in  the  nominal  control  history  required  to  meet  the  final  constraints  while 
minimizing  the  change  in  the  final  time. 

Previous  research  has  successfully  used  this  guidance  law  by  referencing 
the  nominal  control  and  corresponding  change  in  the  control  using  horizontal 
velocity.  This  research  demonstrates  the  guidance  law  using  time  as  the  reference 
variable.  Insuring  that  the  time-to-go  on  both  the  perturbed  and  nominal 
trajectories  are  the  same,  determines  the  change  in  the  control  and  update  of  the 
final  time  that  are  to  be  applied  to  the  perturbed  trajectory. 
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Chapter  1 


Introduction 

Optimal  and  suboptimal  control  theory  are  concerned  with  finding  the 
optimal  control,  and  consequently  the  optimal  trajectory,  for  a  desired  system.  The 
optimal  control  is  found  most  often  using  numerical  optimization.  This  method 
many  times  leads  to  a  less  then  optimal  solution  because  the  model  used  to  describe 
the  system  is  not  accurate.  In  addition,  the  initial  conditions  or  system  parameters 
may  not  be  known  exactly,  or  the  equations  of  motion  may  have  been  simplified. 
As  a  resu't,  when  the  optimal  control  is  applied  to  the  actual  system,  deviations 
from  the  optimal  trajectory  occur,  and  the  optimal  trajectory  is  no  longer  realizable. 
One  solution  to  this  problem  is  to  generate  a  new  optimal  trajectory  from  the  current 
location  to  the  final  constraint  manifold.  However,  because  of  computer 
limitations,  this  solution  is  not  feasible  for  guidance. 

Neighboring  optimal  control  is  another  method  used  to  correct  for 
deviations  from  the  optimal  trajectory.  Assuming  the  deviations  from  the  optimal 
trajectory  caused  by  model  errors  and/or  unknown  perturbations  are  small,  an 
extremal  guidance  law  can  be  developed  to  keep  the  perturbed  trajectory  in  the 
neighborhood  of  the  optimal  or  nominal  trajectory  (see  Ref  1,  for  example). 

This  research  concentrates  on  developing  a  neighboring  extremal  guidance 
law  for  systems  using  piecewise  linear  control  and  demonstrating  the  law  with  time 
as  the  reference  variable.  A  result  of  using  piecewise  linear  control  is  that  the 
control  is  no  longer  a  continuous  function  but  a  set  of  control  parameters.  These 
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control  parameters  arc  the  nodes  used  to  define  the  piecewise  linear  control.  This  is 
also  known  as  suboptimal  control.  The  guidance  law,  using  given  deviations  from 
the  suboptimal  trajectory,  updates  the  control  parameters  using  a  set  of  precomputed 
neighboring  extremal  gains.  The  neighboring  extremal  gains  are  determined  by 
finding  the  change  in  the  control  parameters  that  will  minimize  the  change  in  the 
performance  index  and  still  satisfy  the  final  constraints. 

The  guidance  law  is  demonstrated  using  the  lunar  launch  problem.  The 
lunar  launch  problem  calls  for  placing  a  vehicle  in  a  low  lunar  orbit  in  minimum 
time  while  meeting  conditions  at  orbit  entry.  Perturbations  are  caused  in  the 
trajectory  by  altering  thrust  acceleration  or  lunar  gravity  from  their  nominal  values 
used  in  solving  the  nominal  suboptimal  control  problem.  Previous  work  [2]  has 
successfully  demonstrated  the  guidance  law  using  horizontal  velocity  as  the 
reference  variable  to  determine  the  nominal  control  and  change  in  the  control. 
When  an  attempt  is  made  to  implement  the  guidance  law  with  time  as  the  reference 
variable,  however,  problems  arise.  The  lunar  launch  problem  is  a  free  final  time 
problem  and  the  perturbed  final  time  is  constantly  changing.  As  a  consequence,  it 
is  not  known  exactly  where  the  vehicle  is  located  in  relation  to  the  nominal 
trajectory.  Therefore,  the  nominal  control  and  change  in  the  control  that  must  be 
applied  is  not  known. 

This  research  uses  time  as  the  reference  variable,  and  as  a  result,  the  final 
time  becomes  another  parameter  in  the  control  vector.  The  variation  of  this 
parameter  yields  the  change  in  the  final  time,  which  is  really  the  change  in  the  time- 
to-go  from  a  location  on  the  nominal  trajectory  specified  by  the  reference  time.  The 
change  in  the  time-to-go  is  caused  by  deviations  from  the  nominal  trajectory  at  the 
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location  specified  by  the  reference  time.  The  overall  change  in  the  final  time  is  then 
given  by  the  difference  between  the  reference  time  on  the  nominal  optimal  trajectory 
and  the  current  time  on  the  perturbed  trajectory  such  that  the  time-to-go  on  both 
trajectories  is  the  same. 

Chapter  2  discusses  the  general  optimal  control  problem  and  reformulates 
the  problem  in  terms  of  suboptimal  control.  In  Chapter  3,  the  neighboring 
suboptimal  guidance  law  is  developed.  The  guidance  law  is  then  applied  to  the 
lunar  launch  problem  in  Chapter  4,  with  the  results  shown  in  Chapter  5.  General 
conclusions  and  recommendations  are  discussed  in  Chapter  6. 


Chapter  2 


Optimal  Control 

2.1  Optimal  Control  Problem 

A  dynamical  system  is  governed  by  a  set  of  differential  equations 

x  =  f(t,x,u )  (2.1) 

where  x  is  a  n-vector  of  the  state  variables  and  u  is  a  m-vector  of  the  control.  The 
goal  is  to  guide  the  system  from  some  specified  initial  conditions 

*o  ~  *  *o  =  *o.  (2.2) 

to  the  prescribed  final  conditions 

yr(tf,xf)  =  0,  (2.3) 

by  using  the  control  u.  Many  controls  will  lead  to  trajectories  that  satisfy  these 
boundary  conditions,  but  only  one  is  optimal  in  terms  of  a  desired  performance 
index.  The  optimal  control,  u,  is  determined  by  minimizing  a  performance  index  in 
the  form 

J  =  <j>(t/,x/).  (2.4) 

The  general  optimal  control  problem  can  now  be  stated  as  follows:  find  the  control 
history  u(t)  that  minimizes  Eq.  (2.4),  subject  to  the  differential  equations  given  by 
Eq.  (2.1)  and  the  initial  conditions  of  Eq.  (2.2),  while  satisfying  the  final 
constraints  given  in  Eq.  (2.3)  [1], 

The  solution  to  this  optimal  control  problem  is  usually  too  complicated  to  be 
solved  in  closed  form.  The  solution  to  the  general  optimal  control  problem  may 
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even  be  too  complicated  to  be  solved  numerically.  Therefore,  the  general  optimal 
control  problem  is  reformulated  as  a  suboptimal  control  problem  so  that  a  numerical 
solution  can  be  found. 

2.2  Suboptimal  Control  Problem 

In  the  suboptimal  control,  or  parameter  optimization,  problem  the  control  is 
expressed  as  a  parameter  vector.  This  parameter  vector  contains  the  required 
control  at  different  points  or  nodes  along  the  trajectory.  The  control  between  these 
nodes  is  found  by  interpolating  the  control  between  appropriate  nodes.  In  addition, 
a  free  final  time  problem  can  be  transformed  to  a  fixed  final  time  problem,  and  the 
final  time  becomes  the  last  parameter. 

For  this  research,  consider  the  class  of  piecewise  linear  functions  as  the  set 

of  all  possible  controls.  The  one  control  is  expressed  as  a  parameter  vector 

uT  (2.5) 

where  k  is  the  number  of  control  points  or  nodes  along  the  trajectory.  The  control 

between  two  control  nodes  is  determined  by  linear  interpolating  between  the  two 

control  nodes.  The  continuous  control  function  of  the  optimal  control  problem  is 

now  a  set  of  piecewise  linear  functions. 

Free  final  time  problems  are  transformed  to  fixed  final  time  problems  by 

normalizing  the  variable  of  integration  so  that  the  integration  limits  are  0  to  1.  To 

normalize  the  variable  of  integration,  the  time  is  divided  by  the  final  time,  and  the 

new  variable  of  integration,  t,  is  defined  as 

t 

t  =  — . 


(2.6) 
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The  system  dynamics  can  now  expressed  in  terms  of  the  new  variable  of  integration 
as 

X'  =  ^T  =  t/f(tf'c>x’u)  =  g(t,t/,x,u),  (2.7) 

uT 

where 

T0  =  T0^,  X0  =  ,  Tj  =  1.  (2-8) 

In  this  fixed  final  time  suboptimal  control  problem,  the  states  are  a  function  of  the 
parameterized  control,  u,  and  the  free  final  time.  Thus,  a  different  parameter  vector 
is  defined  which  includes  not  only  the  k  control  nodes  but  also  the  final  time.  This 

new  parameter  vector,  a,  can  be  written  as 

aT  =  [iij,  iij,*",  uk,  tf].  (2.9) 


With  the  new  parameter  vector,  the  object  of  the  fixed  final  time  suboptimal 
control  problem  is  to  find  the  piecewise  continuous  control,  where  xf(a)  is 

obtained  by  integrating  the  system  dynamics 


,  dx  ,  . 

*  =  —  =  g(x,x,a) 
ax 

(2.10) 

subject  to  the  boundary  conditions 

To  =  xo,  >  *0  =  *0.  >  T/  =  1, 

(2.11) 

that  minimizes  the  performance  index 

J  =  <t>(ak+i,x,(aj). 

(2.12) 

and  subject  to  the  boundary  constraints 

v(ak+i,xf(a))  =  0. 

(2.13) 

The  optimal  control  vector,  a,  can  be  obtained  using  a  nonlinear  programming  code 
such  as  VF02AD. 


Chapter  3 


Neighboring  Suboptimal  Control 

The  solution  to  the  suboptimal  control  problem  yields  a  control  history  and 
trajectory  which  minimize  the  desired  performance  index.  For  the  suboptimal 
control  history  to  guide  the  vehicle  to  the  final  boundary  constraints  and  minimize 
the  desired  performance  index,  the  vehicle  must  remain  on  the  nominal  trajectory. 
This  assumption  is  unrealistic,  however,  because  in  order  for  the  vehicle  to  remain 
on  the  nominal  trajectory  the  true  physical  model  and  the  model  used  to  generate 
the  suboptimal  control  would  have  to  be  exactly  the  same.  In  reality,  the  vehicle 
will  not  remain  on  the  nominal  trajectory  because  of  perturbations  in  the  initial 
conditions  and/or  model  errors.  As  a  consequence,  a  guidance  law  is  required  that 
operates  the  vehicle  in  the  neighborhood  of  the  suboptimal  trajectory. 

Neighboring  suboptimal  control  is  used  to  develop  a  guidance  law  for 
operating  a  vehicle  in  the  neighborhood  of  the  suboptimal  path.  Given  that  a 
perturbation  away  from  the  desired  states  occurs,  the  goal  of  neighboring 
suboptimal  control  is  to  determine  the  change  in  the  nominal  control  required  so 
that  the  overall  change  in  the  performance  index  is  minimized  and  tne  boundary 
conditions  are  satisfied. 

The  goal  of  developing  a  neighboring  suboptimal  control  law  is  to  find  a 
relationship  between  deviations  from  the  nominal  trajectory  and  changes  required 
in  the  nominal  control  so  that  the  change  in  the  desired  performance  index  is 
minimized  and  the  boundary  conditions  are  satisfied.  The  derivation  of  this 
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control  law  was  originally  developed  by  Hull  and  Helfrich  [3].  A  summary  of  the 
derivation  is  included  here  for  completeness. 

As  shown  in  Equations  (2.12)  and  (2.13),  the  value  of  the  performance 
index  and  the  constraints  are  a  function  of  the  parameter  vector  a.  Therefore,  the 
augmented  performance  index  for  the  suboptimal  control  problem  is  expressed  as 
S'  =  <fi(ak+1,x,(a))+  vTv(aM  ,x,(a))  =  G(at+1,X/(a),v)  =  G(a,v),  (3.1) 

where  0  is  the  performance  index,  y/  is  the  vector  of  final  constraints  and  v  is 
the  vector  of  Lagrange  multipliers.  Given  a  perturbation  in  the  control  vector,  a, 
the  change  in  the  performance  index  is 

A/'  =  G(a  +  8a,v)-  G{a,  v).  (3.2) 

Since  the  first  variation  of  the  performance  index  vanishes  on  the  optimal  path, 
the  change  in  the  performance  index  due  to  a  perturbation  in  the  control  vector  to 

second  order  is  approximated  as 

A  J'  =  ±8aTGJa,v)5a.  (3.3) 

The  second  derivative  matrix  GM  is  calculated  numerically  for  the  optimal  path. 
In  Eq.  (3.3),  the  5d  s  are  not  independent  but  constrained  by  the  requirement  that 
the  final  constraints  be  satisfied.  Therefore,  a  relationship  is  needed  to  correlate  a 
perturbation  in  the  control,  5a,  to  a  deviation  from  the  optimal  path,  <5x0,  such 
that  5y/=  VXf8xf  =0. 

The  system  dynamics  given  by  Eq.  (2.10)  relate  the  state  and  the  control 
through  a  set  of  differential  equations.  To  relate  the  variation  8a  to  the 
perturbation  8x0,  the  variation  of  Eq.  (2.10)  is  taken  as 

8x'  =  gx8x  +  ga5a. 


(3.4) 
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with  boundary  conditions 


II 

o'* 

h 

II 

►— * 

(3.5) 

#• 

II 

* 

«• 

(3.6) 

yrX/8xf=0. 

(3.7) 

To  solve  Eq.  (3.4),  a  solution  is  assumed  to  be  of  the  form 

&  =  0(T,l)&y  +'P(t,I)& 

(3.8) 

with  boundary  conditions 

O(U)  =  I 

(3.9) 

mi)=o. 

(3.10) 

Substituting  the  derivative  of  Eq.  (3.8)  into  the  left  side  of  Eq.  (3.4),  and  Eq.  (3.8) 

into  the  right  side  of  Eq.  (3.4)  results  in 

(O'  -  gx<S>)8xf  +  (¥'  -gx^-ga)8a  =  0.  (3. 1 1 ) 

To  guarantee  that  Eq.  (3.11)  is  satisfied,  the  coefficients  of  Sxf  and  8a  are  chosen 
to  be  zero.  The  outcome  of  setting  the  coefficients  equal  to  zero  is  two 
differential  equations  defining  the  transition  matrices  O  and  T*.  The  equations 
are 

0'  =  g,0  (3.12) 

and 

(313) 

with  boundary  conditions 

*/  =  I 
^=0. 


(3.14) 

(3.15) 
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By  using  Eqs.  (3.12)  and  (3.13),  O'  and  can  be  integrated  backwards  along  the 
nominal  trajectory  from  rf  to  find  <t>(r)  and  *F(t).  Knowledge  of  the  history  of 

the  transition  matrices  allows  Eq.  (3.8)  to  be  used  at  any  time. 

Recall  the  variation  of  the  system  dynamics  ,  given  by  Eq.  (3.4),  is  subject 
to  the  constraint  of  Eq.  (3.7).  Multiplying  Eq.  (3.8)  by  O"1  and  solving  for  Sxf 

yields 

Sxf  =  -d +  Q-'Sx.  (3. 16) 

Applying  Eq.  (3.16)  at  f0  and  substituting  this  result  into  Eq.  (3.7)  leads  to 

V^oX&^O,  (3.17) 

which  is  a  constraint  on  5a  imposed,  for  a  given  5x0,  as  a  result  of  satisfying  the 
final  conditions. 

A  new  optimization  problem  is  now  formed  to  minimize  the  change  in  the 
performance  index  given  by  Eq.  (3.3)  with  respect  to  5a,  subject  to  the  constraint 
of  Eq.  (3.17).  The  new  augmented  performance  index  is  written  as 

f  =  jj5aTGaa5a  +  5vT[yfXi<S>-j8x0  -  y^d^^Sa).  (3.18) 

By  minimizing  Eq.  (3.18),  the  change  in  the  performance  index  is  minimized  due 
to  perturbations  in  the  control  and  subject  to  the  boundary  constraint  that  relates 
the  control  and  the  states. 

Taking  the  first  variation  of  Eq.  (3. 18)  results  in 

Sy  =  (SaTGm  +  SvTy,I®~0"¥0)S(Sa).  (3.19) 

To  find  a  minimum,  the  first  variation  must  vanish.  This  implies  that 

&rGM  +  5vVX/d>-X  =  0 


(3.20) 
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needs  to  be  satisfied  for  a  minimum  to  exist.  Transposing  Eq.  (3.20)  and 
assuming  exists  leads  to 

8a  =  -G2'Vl<t>oTvlfSv.  (3.21) 

Substituting  Eq.  (3.21)  into  Eq.  (3.17)  yields 

-  y/Xf  O-I'F0G-1'Fj<D"r  y/TXf  Sv  =  y/X/  (3.22) 

which  can  be  solved  for  Sv  as 

Sv^-(¥,0l"V,G^l€';TwTx,j\.,^'Sxo-  (3.23) 

Inserting  Sv  into  the  right  side  of  Eq.  (3.21)  leads  to  the  neighboring  suboptimal 
control  law 

8a  =  K08x0,  (3.24) 

where 

=  Ci'W  (  Vl>  OiXGi'W  < )" V,, 4>;‘  (3.25) 

is  the  gain  matrix  at  r0. 

The  concept  of  controllability  assumes  that  to  first  order  a  trajectory  exists 
from  a  perturbed  location,  Sx0,  at  a  given  time,  T,  to  the  final  constraint  manifold, 

or  Sy/  ~  0.  In  other  words,  a  solution  to  Eq.  (3.4)  is  desired  which  satisfies  the 
boundary  conditions  given  by  Eqs.  (3.5-3.7).  Recall  that  to  solve  for  the 
suboptimal  control  gains  given  in  Eq.  (3.25),  a  solution  to  Eq.  (3.4)  is  assumed  in 
the  form  of  Eq.  (3.8).  Therefore,  the  system  is  assumed  to  be  controllable.  This 
implies  there  exists  at  least  one  control  history  that  will  lead  from  the  current 
perturbed  location  to  the  final  constraint  manifold.  The  neighboring  suboptimal 
control  law  determines  the  one  control  history  that  will  minimize  the  change  in 
the  performance  index  while  satisfying  the  final  constraints.  Thus,  the  control  law 
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given  in  Eq.  (3.24)  uses  small  deviations  in  the  states,  Sx0,  at  a  given  time,  r,  to 

determine  small  changes  required  in  the  parameters,  8a.  These  small  changes  are 

then  added  to  the  optimal  values  of  the  control  parameters,  found  by  solving  the 
suboptimal  control  problem  from  r0  to  rf,  to  form  a  new  control  history.  This 

new  control  history  leads  to  the  final  constraints  while  minimizing  the  change  in 
the  performance  index. 

This  method  of  control  can  be  extended  so  that  any  time  Ti<T/, 
i  =  may  be  considered  t0  or  the  "initial"  time  and  the  control  history 

updated.  Because  <J>  and  'F  are  functions  of  r  only,  the  gain  K0  can  be  computed 
and  stored  at  each  node  point  along  the  nominal  optimal  trajectory.  When  r0 
occurs  between  node  points,  the  gain  is  linearly  interpolated.  Thus,  at  each 
sample  time  the  current  states  of  the  vehicle  are  compared  with  the  states  along 
the  nominal  optimal  path,  the  appropriate  gains  chosen,  and  the  change  in  the 
control  parameters  determined. 


Chapter  4 


Application  of  the  Control  Law 


4.1  Lunar  Launch  Problem 

The  neighboring  suboptimal  control  law  is  demonstrated  using  a  lunar 
launch  problem.  The  lunar  launch  problem  involves  placing  a  vehicle  into  lunar 
orbit  in  minimum  time  (refer  to  Fig.  4.1).  In  addition,  constraints  are  placed  on 
the  vehicle's  position  and  velocity  at  orbit  entry.  The  control  used  to  place  the 
vehicle  in  orbit  is  the  angle  of  thrust  versus  time,  6(t).  The  lunar  launch  problem 
is  a  free  final  time  problem  whose  solution  provides  a  control  history  which 
minimizes  the  performance  index 

J  =  tf.  (4.1) 

The  differential  equations  governing  the  trajectory  of  the  vehicle  are  given 
by 

x  =  u  (4.2) 

y  =  v  (4.3) 

ii=  a  cos  0  (4.4) 

v=  asin0-g,  (4.5) 

where  a  is  the  thrust  acceleration  and  g  is  the  acceleration  due  to  gravity.  The 
prescribed  initial  conditions  are 

to  =  0,  *0  =0,  y0  =0,  «o  =0,  v0  =0,  (4.6) 

and  the  final  conditions,  or  the  conditions  at  orbit  entry  are 
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Figure  4. 1  The  lunar  launch  problem 


yf  =  50,000  ft 

(4.7) 

uf  =  5,444  ft  /sec 

(4.8) 

v/=0  ft  /sec. 

(4.9) 

The  quantities  u  and  v  denote  the  horizontal  and  vertical  velocities,  respectively. 

To  simplify  the  problem,  several  assumptions  have  been  made.  One 
assumption,  already  implied  by  the  equations  of  motion,  is  that  all  motion  and 
thrust  act  in  the  x-y  plane  only.  In  addition,  the  vehicle’s  trajectory  is  assumed  to 
be  over  a  flat  moon,  and  therefore  the  direction  of  lunar  gravity  remains  constant. 
Also,  the  magnitude  of  lunar  gravity  is  assumed  to  remain  constant  throughout  the 
trajectory.  Finally,  the  ratio  of  vehicle  mass  and  thrust  is  assumed  constant 
throughout  the  trajectory.  This  implies  that  thrust  acceleration,  a,  is  constant. 
Based  on  previous  research  on  the  lunar  launch  problem  [31,  values  of 
g  =  5.32 ft  / sec2  and  a  =  20. Sft  /  sec2  are  assumed. 
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The  current  formulation  of  the  lunar  launch  problem  is  as  a  free  final  time 
problem  and  calls  for  a  continuous  control  history  solution.  By  using  the  ideas  of 
suboptimal  control,  T  is  introduced  as  the  new  variable  of  integration  and  is 
normalized  with  respect  to  the  final  time.  The  normalized  differential  equations 
become 


,  dx  dx  dt  t  „  \ 

x  =—=——= g(T,x,  e,tf\. 
dx  dt  dx  '  ' 

(4.10) 

The  equations  of  motion  for  the  vehicle  are  now 

x  =  tfu 

(4.11) 

> 

♦Jr- 

ii 

(4.12) 

ii=  rfacos0 

(4.13) 

v=  rf(a  sin  0-g). 

(4.14) 

The  initial  conditions  remain 

xT0  =[0  0  0  0] 

(4.15) 

and  the  limits  of  integration  become 

it 

o 

4"* 

II 

M 

(4.16) 

To  find  the  suboptimal  control  history,  nine  control  nodes  are  evenly 
spaced  throughout  the  optimal  trajectory.  The  control  applied  at  each  node,  0i , 
and  the  optimal  final  time,  tf,  are  used  to  define  the  suboptimal  control,  or 

parameter,  vector,  a.  The  optimal  control  vector,  a ,  is  defined  as 

aT=[e{,  02, -,  09,  tf\  (4.17) 

The  final  conditions,  normalized  for  numerical  reasons  to  have  their 
magnitudes  of  the  same  order,  are  placed  in  a  constraint  vector,  y/ ,  where 
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-2L-.-1 

50000 


V  = 


1 

5444 


=  0. 


5444 


(4.18) 


The  optimal  control  vector,  a ,  is  computed  using  VF02AD,  a  nonlinear 
programming  code  [4].  The  VF02AD  program  is  an  extension  of  a  variable 
metric  method  for  unconstrained  optimization  applied  to  constrained  optimization. 
The  program  uses  a  recursive  quadratic  programming  variable  metric  method.  A 
suboptimal  solution  is  searched  for  until  a  preset  convergence  criteria,  in  this  case 
10”9  sec,  is  reached.  Convergence  occurs  when  the  change  in  the  performance 
index  between  iterations  is  less  then  the  preset  convergence  criteria. 

4.2  Calculation  of  the  Neighboring  Suboptimal  Control  Gains 


Recall  that  the  gain  K0  is  a  function  of  x  and  is  multiplied  by  the 

deviation  from  the  nominal  optimal  trajectory  to  determine  the  required  change  in 
the  nominal  control.  Since  K0  is  defined  as 

<419> 

to  calculate  K0  the  differential  equations  for  O  and  'P  must  be  solved  and  the 
expressions  for  and  y/Xf  evaluated. 

To  solve  for  O  and  H*,  their  corresponding  differential  equations  must  be 
integrated  backwards  from  rf.  The  differential  equation  governing  <J>  is  given  by 

~  =  4>'  =  gx<t>  (4.20) 

ax 


where 
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&f=I.  (4.21) 

The  expression  gx  is  evaluated  as  the  partial  of  the  system  dynamics,  given  in 


Eqs.  (4.11-4.14),  with  respect  to  the  states  and  is  given  by 

'0  0  tf  O' 

dg  0  0  0  ty 

gl  =  dx=  0  0  0  0' 

0  0  0  0 


(4.22) 


The  differential  equation  governing  4*  is 

'T  =  gx'¥  +  ga  (4.23) 

where 

¥,=0.  (4.24) 

In  Eq.  (4.23),  the  expression  for  gx  is  given  by  Eq.  (4.22).  The  expression  ga  in 


Eq.  (4.23)  represents  the  partial  derivative  of  the  system  dynamics  with  respect  to 

the  parameter  vector  a.  Because  a  is  a  vector  and  is  expressed  as 

aT  =  [dv  02,~  ,  09,  tf],  (4.25) 

the  chain  rule  must  be  applied  to  form  gB .  Since  there  are  10  parameters,  ga  is 


written  as 


-  dg  da 
ga  da  da: 


1  Si  <  10. 


(4.26) 


The  partial  derivatives  determine  the  impact  a  change  in  one  parameter  has  on 
another.  The  partial  derivative  relating  the  change  in  the  system  dynamics  to  the 


control  vector,  a,  can  be  expressed  by 


dg  _  dg  Qr  dg 
da  d6  dt, 


(4.27) 


The  expressions  are  evaluated  as 
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dg__ 

de 


and 


dg  _ 

dt. 


dx 

de 

dy 

de 

du 

de 

dv 

del 


dx 

dt, 

dy 


0 

0 

-asin&t, 
a  cos  6tf 


(4.28) 


du 

*7 

!k. 


u 

v 

a  cos  0 
a  sin#  -g 


(4.29) 


To  calculate  the  partial  derivative  of  the  parameter  vector,  a ,  with  respect 
to  each  element  of  the  parameter  vector,  at,  recall  the  optimal  control,  given  in  the 

first  9  elements  of  a ,  is  linearly  interpolated  between  adjoining  nodes  as 

0(f)  =  0f._,  +  — — — (t-  t,_,),  for  T(-_]  <  T <  T,  (4.30) 


T-  -  T  , 


and 


0(r)  =  0,  +  — — — (t-  t(),  for  r,  <  t <  t, 


/+!• 


(4.31) 


Therefore,  the  partial  derivatives  of  Eqs.  (4.30)  and  (4.31)  with  respect  to  0,, 
respectively,  are 


00 

00. 


7,-1 


T;  ~  T; 


for  T,_,  <  T  <  T, 


(4.32) 


-1 


and 
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90  T—  T, 
~  fi+I  -  T, 


for  r,  <  t  <  t,+1  . 


(4.33) 


Both  Eqs.  (4.32)  and  (4.33)  are  valid  for  all  interior  node  points.  However,  each 
equation  breaks  down  at  one  of  the  boundary  nodes  because  either  rt_,or  ti+1 


does  not  exist.  Therefore,  only  Eq.  (4.33)  is  used  over  the  first  interval  and  only 
Eq.  (4.32)  is  used  over  the  last  interval.  The  last  element  of  a  is  tf,  and  the  partial 
of  tf  with  respect  to  tf  is  one.  Over  all  other  intervals,  the  partial  derivative  is 


zero. 


The  matrix  ga  is  expressed,  using  Eqs.  (4.28),  (4.29),  (4.32)  and  (4.33),  as 

0  •••  0  u 

0-0  v 


ga  = 


.  Q  dd  .  _  dd  a 

-asm  Or,——  —  -asm  Or,— —  acosO 

f  dd ,  1  de9 


acosOt, 


90 


a  cos  Or, 


90 


a  sinO-  g 


(4.34) 


f  90 ,  f  909 

The  derivative,  GB,  represents  the  second  derivative  effect  of  a  change  in 
the  parameter  vector  on  the  performance  index  and  can  be  determined 
numerically.  Central  differences  are  employed,  and  each  element  of  a  is 
perturbed.  The  effect  of  this  perturbation  on  the  performance  index  is  computed. 

The  partial  derivative  of  the  constraints  with  respect  to  the  state  vector  is 
represented  by  y/x  .  This  forms  a  constant  matrix  given  by 


0 


1 


50,000 
0  0 


0 

1 


0 


0 


5,444 

0 


0 

0 

1 


5,444 


(4.35) 
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By  using  Eq.  (4.35),  the  values  computed  for  G m,  and  the  solution  to  the 
differential  equations  for  d>  and  4*,  the  general  gain  matrix  given  in  Eq.  (4.19) 
can  be  solved  at  any  control  node.  These  gains,  multiplied  by  the  deviation  away 
from  the  nominal  trajectory  at  the  corresponding  node  point,  yield  the  required 
change  in  the  control.  Because  the  first  column  in  Eq.  (4.35)  is  all  zeros,  the  first 
column  of  the  gain  matrix  will  also  be  all  zeros.  This  is  not  surprising  because  the 
downrange  is  not  constrained,  and  therefore  any  deviation  in  the  downrange  from 
the  nominal  trajectory  will  not  affect  the  control.  This  contrasts  all  of  the  other 
states,  where  any  deviation  from  the  nominal  trajectory  will  require  a  change  in 
the  control  to  insure  the  final  constraints  are  satisfied.  As  a  result,  only  the  last 
three  columns  of  the  gain  matrix  will  be  stored. 

Eight  gain  matrices  are  computed  and  correspond  to  the  first  eight  control 
nodes.  The  last  control  node  is  at  the  final  constraint,  and  therefore  a  gain  matrix 
can  not  be  computed  because  there  are  not  enough  nodes  to  satisfy  the  boundary 
conditions.  If  the  gain  matrix  associated  with  the  ninth  node  could  be  computed, 
however,  it  would  call  for  an  infinite  control  change  to  correct  for  any  deviations 
from  the  final  constraints.  As  a  result,  the  gains  from  the  previous  two  nodes  will 
be  used  to  extrapolate  the  change  in  the  control  over  the  last  interval.  The  gain 
matrix,  Ki ,  contains  information  to  change  the  value  of  the  control  at  node  i  and 

all  nodes  after  node  i,  given  a  deviation  from  the  nominal  trajectory  occurred  at 
node  i.  The  control  law  will  be  applied  throughout  the  trajectory,  however,  and 
only  the  gains  associated  with  node  i  and  the  node  immediately  following  node  i 
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are  needed  for  applying  the  control  law.  Therefore,  only  these  gains  will  be  stored 
from  the  gain  matrix  Ki . 


4.3  Interpolating  the  Change  in  the  Control 


The  eight  gain  matrices  provide  the  required  changes  in  the  control  nodes 
from  node  i  to  the  final  constraint  manifold  for  a  deviation  from  the  nominal 
trajectory  occurring  at  node  i.  The  change  in  the  control  nodes  are  then  added  to 
the  nominal  controls  to  provide  a  new  control  history.  If  the  deviation  from  the 
nominal  trajectory  occurs  between  nodes,  the  gains  from  two  nodes  are  used  and 
the  change  in  the  control  linearly  interpolated  between  the  two  nodes  to  form  a 
new  control  history.  Previous  research  by  Nowack  [2]  has  successfully  used  a 
sample  hold  technique  for  interpolating  the  change  in  the  control  between  nodes. 
The  nominal  control,  however,  is  linearly  interpolated  throughout  the  sample  step. 
Therefore,  a  new  method  of  interpolating  the  change  in  the  control  is  implemented 
which  will  linearly  interpolate  the  change  in  the  control  throughout  the  sample 
step. 


Assuming  a  deviation  from  the  nominal  trajectory  occurred  at  node  i,  the 
change  in  the  control  nodes  after  node  i  are  given  by 


8ei.i  =  Ki.yfyi +  Kj.M  + 


(4.36) 


where 

i  =  current  node  where  state  deviations  occur 
j  =  node  where  the  change  in  the  control  is  applied 
8yt  =  deviation  in  y  from  nominal  trajectory 
8Uj  =  deviation  in  u  from  nominal  trajectory 
<5v;  =  deviation  in  v  from  nominal  trajectory 
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and  the  gain  matrix  is  partitioned  as 


k:  .  v 

Ic  k 

»v;  1  fv;  | 

K.l.y 

1,1,1*  v 

k  k 

*1,2,1*  *i,  2,v 

"<Sy. 

5m,. 

* 

•  * 

A. 

Kj.y 

k-  ■  k- 

*i,/^*  A'*,/,v_ 

The  expression  ;  is  only  valid  when  i  <  j ,  because  it  otherwise  implies 
changes  to  control  nodes  that  have  already  been  passed.  Also,  notice  that  all 
terms  involving  x  have  been  dropped  because  any  deviation  in  this  state  would  not 
affect  the  control.  Realize  also  that  Eq.  (4.37)  uses  deviations  in  three  states  to 
determine  the  change  in  the  control.  This  is  in  contrast  to  Nowack  [2]  in  which 
horizontal  velocity,  u,  is  used  as  the  reference  variable.  No  deviation,  5m,,  ever 

exists,  and  therefore  only  deviations  in  y  and  v  are  used  to  determine  the  change  in 
the  control. 

The  change  in  the  control  nodes,  given  in  Eq.  (4.36),  are  added  to  the 
nominal  controls,  00  ,  at  each  of  the  remaining  nodes.  The  new  control  history  is 

formed  by  linearly  interpolating  between  the  new  control  nodes,  ,  where 

ej  =  0Oj+seu.  (4.38) 

If  the  deviation  from  the  nominal  trajectory  occurs  between  nodes, 
however,  no  explicit  gain  matrix  exists  and  the  expression  80^  cannot  be 

computed  directly.  Therefore,  the  gain  matrices  corresponding  to  the  control 
nodes  before  and  after  the  deviation  must  be  used  to  interpolate  the  appropriate 
change  in  the  control.  If  the  deviation  from  the  nominal  trajectory  had  occurred  at 
node  i,  a  new  control  history  could  be  determined  from  node  i  to  the  final 
constraint  manifold  using  50,  y.  Similarly,  if  the  deviation  had  occurred  at  node 

i+1,  a  different  control  history  from  node  i+l  to  the  final  constraint  manifold 
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Figure  4.2  Change  in  the  control 


could  be  determined  using  S0M  j.  It  seems  reasonable,  therefore,  that  for  a 

deviation  occurring  between  two  nodes  there  exists  another  control  history  that 
lies  somewhere  between  the  two  control  histories.  This  new  control  history  is 
likewise  determined  by  using  a  change  in  the  control  that  lies  somewhere  in 
between  SQ^  and  S0MJ. 

The  change  in  the  control  for  a  deviation  occurring  at  time  t,  between 
nodes  1  and  2  (see  Fig  4.2),  can  be  determined  by  first  assuming  that  the  deviation 
from  the  nominal  trajectory  occurred  at  the  previous  node,  node  1,  and  by  then 
computing  the  corresponding  change  in  the  control,  5du,  where 

Mu=ku.,fy  +kujM  +kui,M-  <4-39> 

Now,  assume  the  deviation  occurred  at  the  next  node,  node  2,  and  compute  the 
change  in  the  control,  S022,  from  node  2  to  the  final  constraint  manifold  using 

802,2  =  *2,2,  yfy  +  k2.2,uSu2  +  *2.2.v5v2  *  C4-4*)) 
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To  determine  the  appropriate  change  in  the  control  between  nodes,  the 
sample  hold  technique  used  by  Nowack  [2]  calls  for  linear  interpolating  between 
80u  and  S622  to  find  86 A.  This  quantity,  given  by 

(4.4i) 

h  h 

is  then  held  constant  throughout  the  sample  step  and  added  to  the  nominal  control. 
The  nominal  control,  however,  is  linearly  interpolated  throughout  the  sample  step. 
Therefore,  to  also  linearly  interpolate  the  change  in  the  control  throughout  the 
sample  step,  the  quantity  86 B  is  calculated.  This  quantity  is  found  by  linearly 
interpolating  the  change  in  the  control  between  the  control  nodes  which 
immediately  follow  the  two  nodes  where  the  deviation  is  assumed  to  have 
occurred.  In  this  case  the  change  in  the  control  required  at  node  2,  due  to  a 
deviation  occurring  at  node  1,  is  given  by  86h2.  The  change  in  the  control 
required  at  node  3,  due  to  a  deviation  occurring  at  node  2,  is  given  by  8023.  This 
quantity  86 B  (see  Fig.  4.2)  is  found  by  linearly  interpolating  between  80l  2  and 
S02  3  and  is  given  by 

se,  =  66,  ^  (4.42) 

h  h 

The  two  quantities  <50Aand  86 B  are  now  used  to  define  the  change  in  the  control 
required  for  a  deviation  from  the  nominal  trajectory  occurring  at  time  t  between 
nodes  1  and  2.  During  the  sample  step,  the  change  in  the  control  is  linearly 
interpolated  between  86 A  and  86 B  as 

86  =  86 A  +  (/-/,)  (4.43) 

and  added  to  the  nominal  control  to  form  the  new  control. 
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This  technique  could  be  used  to  determine  the  entire  new  control  history 
from  the  point  of  the  deviation  to  the  final  constraint  manifold,  but  this  is  not 
required.  At  each  sample  time,  the  previous  control  history  would  be  discarded 
and  a  new  control  history  computed.  Therefore,  the  change  in  the  control  history 
is  computed  only  between  the  two  nodes  needed  to  linear  interpolate  the  change  in 
the  control. 

The  last  control  node  is  located  at  the  final  constraint  manifold.  There 
does  not  exist  a  gain  matrix  at  this  node.  As  a  result,  the  gain  matrices  from  the 
previous  two  nodes  are  extrapolated  over  the  last  interval.  The  ideas  and  concepts 
previously  discussed  remain  valid,  but  the  equations  are  slightly  altered  to  linearly 
extrapolate,  instead  of  interpolate,  over  the  last  interval. 

4.4  Calculation  of  the  Reference  Time 

AH  the  expressions  are  now  known  to  solve  the  gain  matrix  at  each  node. 
The  interpolation  of  the  gain  matrices  between  nodes  to  determine  the  change  in 
the  control  given  a  specific  reference  time,  t0,  is  also  understood.  It  is  not  known, 
however,  at  what  specific  time,  t0,  the  gains  are  to  be  interpolated.  Recall  that  the 

gains,  as  well  as  the  nominal  control,  are  evaluated  with  respect  to  T  and  that  r  is 
related  to  both  a  time  on  the  nominal  trajectory  and  the  optimal  final  time,  tf  . 

Assuming  the  initial  conditions  are  the  same  at  the  beginning  of  the  trajectory,  the 
final  time  of  the  perturbed  trajectory  is  assumed  the  same  as  the  optimal  final  time 
of  the  nominal  trajectory.  But,  as  soon  as  some  deviation  from  the  nominal 
trajectory  occurs,  the  optimal  final  time  is  no  longer  realizable.  If  the  current  time 
on  the  perturbed  trajectory  is  used  to  determine  the  appropriate  gain  to  apply,  the 
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obvious  problem  arises  when  the  final  time  of  the  perturbed  trajectory  is  greater 

than  the  final  time  of  the  nominal  trajectory  [5,6].  The  gains  and  the  nominal 
control  are  known  only  until  r  =  1,  which  corresponds  to  a  time  equal  to  tf  . 

Therefore,  to  insure  that  both  gains  and  a  nominal  control  are  available  throughout 
the  trajectory,  a  method  which  incorporates  time-to-go  is  used  to  reference  both 
the  nominal  control  and  the  gains  associated  with  the  change  in  the  nominal 
control. 

The  lunar  launch  problem  is  a  free  final  time  problem  whose  nominal 
solution  yields  an  optimal  final  time,  control  history,  and  trajectory  to  the  final 
constraint  manifold.  The  controls  are  applied  by  determining  the  location  of  the 
vehicle  on  the  nominal  trajectory  with  respect  to  T.  If  the  vehicle  remains  on  the 
nominal  trajectory,  the  final  time  does  not  change,  and  the  current  time  is  used  to 
reference  the  control.  When  the  vehicle  deviates  from  the  nominal  trajectory, 
however,  the  final  time  changes  and  the  location  of  the  vehicle  with  respect  to  the 
nominal  trajectory  is  unknown.  Thus,  it  is  not  known  what  gains  should  be  used 
to  determine  the  change  in  the  control,  or  even  what  nominal  control  should  be 
applied.  Therefore,  using  the  ideas  of  Kelley  [5]  and  Speyer  and  Bryson  [6],  a 
location  on  the  nominal  trajectory  is  found  such  that  the  time-tc-go  both  on  the 
nominal  trajectory  and  the  perturbed  trajectory  is  the  same.  This  provides  a 
location  of  the  vehicle  with  respect  to  the  nominal  trajectory  and  a  reference  time 
to  determine  the  nominal  control  and  the  change  in  the  control.  Having  both  the 
nominal  and  perturbed  trajectories  always  referenced  at  the  same  time-to-go 
guarantees  that  a  nominal  control  and  change  in  the  control  will  be  available 
throughout  the  entire  trajectory,  regardless  of  the  perturbed  final  time. 
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Figure  4.3  Relating  time-to-go 

The  change  in  the  final  time  caused  by  a  deviation  from  the  nominal 
trajectory  is  expressed  as 

= '-v  <444> 

where  is  the  perturbed  final  time,  tf^  is  the  optimal  final  time  on  the  nominal 

trajectory,  t  is  the  actual  time  on  the  perturbed  trajectory  and  t0  is  a  reference  time 

on  the  nominal  trajectory  (see  Fig  4.3).  Equating  time-to-go  on  the  perturbed  and 

nominal  trajectories  results  in 

tfr~t  =  tf<m  - 10  +  <5(r^  - 10 ) ,  (4.45 ) 

where  S[tf^  -  f0)  is  the  difference  between  time-to-go  on  the  perturbed 

trajectory  and  the  time-to-go  on  the  optimal  trajectory  at  the  reference  time,  t0. 
Recall  the  last  element  of  the  parameter  vector,  a,  is  defined  as  tf .  The  variation 

of  this  expression  defines  the  change  in  the  final  time,  where  the  final  time  is  the 
time  from  the  reference  time,  t0 ,  to  the  optimal  final  time,  tf  .  The  gains  used  at 
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the  reference  time,  r0,  were  evaluated  by  integrating  the  transition  matrices,  given 
by  Eqs.  (3.12-3.15),  backwards  from  r/,  or  tf ,  to  the  control  node  where  the 

gains  are  stored.  The  control  nodes  are  identified  in  each  integration  as  T0,  but 

since  the  optimal  final  time  is  known,  the  control  nodes  could  also  be  referenced 
by  their  corresponding  reference  time,  t0.  Thus,  Sty  relates  a  change  in  time,  not 

over  the  entire  trajectory,  but  from  the  reference  time,  t0,  to  the  end.  Therefore, 
the  quantity  8tf  can  be  written  as  S(tf^  -  r0j,  or  the  change  in  the  time-to-go 

from  the  nominal  trajectory  to  the  perturbed  trajectory. 

The  goal  is  to  find  a  reference  time,  t0 ,  where  the  time-to-go  on  each 
trajectory  is  the  same.  This  implies  that  Stf  is  equal  to  zero.  If  r0  is  the  reference 

time  where  the  time-to-go  on  both  trajectories  is  the  same,  then  a  time  earlier  than 
t0  will  yield  a  positive  8tf  and  a  time  after  t0  a  negative  8tf  (see  Fig.  4.3).  Thus, 

a  one-dimensional  search  routine  is  used  to  search  for  a  reference  time  where  8tf 

is  zero.  Because  the  routine  searches  for  a  minimum,  the  input  for  the  search  is 
the  reference  time  from  the  previous  step  and  a  search  is  made  over  (&y)  .  The 

output  of  the  one-dimensional  search  is  a  new  reference  time  on  the  nominal 
trajectory.  This  new  reference  time  corresponds  to  the  new  time  on  the  perturbed 
trajectory  and  indicates  where  8tf  is  zero.  This  new  reference  time  is  then  used  to 

determine  the  nominal  control  and  the  change  in  the  control  required  to  satisfy  the 
final  constraints.  The  difference  in  the  actual  time  on  the  perturbed  trajectory,  /, 
and  the  new  reference  time,  r0,  furnishes  a  new  estimate  of  the  perturbed  final 

time  using  Eq.  (4.44). 
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Forcing  the  time-to-go  on  the  perturbed  trajectory  to  be  the  same  as  the 
time-to-go  on  the  nominal  path  assumes  that  there  are  no  further  perturbations  in 
the  system.  If  the  previous  perturbations  in  the  system  are  due  to  errors  in  the 
model,  they  will  not  be  accounted  for  in  the  new  control  history.  The  new  control 
history  is  based  on  the  same  model  used  to  generate  the  optimal  control  vector. 
Thus,  the  actual  perturbed  final  time  is  not  known  until  it  is  actually  reached. 
Throughout  the  trajectory,  the  process  continues  to  update  the  perturbed  final  time 
and  moves  it  closer  to  the  actual  final  time. 

4.5  Simulation 

A  simulation  using  the  lunar  launch  problem  is  made  to  test  the  guidance 
law.  The  differential  equations  of  motion  for  the  vehicle  are  integrated  using  a 
fourth-order  Runge-Kutta  integrator  and  a  4  second  sample  step.  Deviations  from 
the  nominal  trajectory  are  introduced  by  perturbing  the  values  of  thrust 
acceleration  or  lunar  gravity  ±5%  from  their  nominal  values.  The  nominal 
trajectory  is  integrated  to  the  appropriate  reference  time  and  deviations  in  the 
states  are  measured.  The  simulations  are  accomplished  closed  loop  using  both  the 
linear  interpolated  control  and  change  in  the  control.  The  results  are  compared 
with  simulations  obtained  by  Nowack  [2]  in  which  the  reference  variable  was 
changed  to  horizontal  velocity. 


Chapter  5 


Results 

5.1  Nominal  Suboptimal  Trajectory 

The  nine  node  suboptimal  trajectory  and  corresponding  final  time  are  first 

calculated  using  the  nonlinear  programming  code  VF02AD.  The  result  of  this 

yields  the  following  control  history: 

■AH  [0. 45427200240]  [26.02786848706' 

02  0.40895291023  23.43127577463 

03  0.36254811487  20.77247685115 

04  0.31389511079  17.98486505714 

05  =  0.26421940074  rad  =  15.13865652683  deg. 

06  0.21307434261  12.20826055367  (5  1) 

07  0.16035307712  9.18755455052 

0*  0.10711091552  6.13700339899 

09J  [0.05217078 106  J  [  2.98916556845 

The  vehicle  follows  the  nominal  trajectory  using  the  control  history  given  in  Eq. 

(5.1)  and  reaches  orbit  with  the  final  conditions: 


tf  =  272. 70605  sec 

J  opt 

(5.2) 

Xj  =  729,138.05373  ft 

(5.3) 

yf  =  50,000.00000  ft 

(5.4) 

«/=  5,444.00000  ft/sec 

(5.5) 

vf  =0.00000  ft/sec. 

(5.6) 

30 


31 


5.2  Neighboring  Suboptimal  Control  Gains 

The  gain  calculations  are  made  using  the  nominal  control  history  given  in 
Eq.  (5.1).  The  gain  program  first  integrates  the  equations  of  motion  to  the  final 
constraint  manifold  and  then  integrates  backwards  to  each  of  the  control  nodes. 
Only  certain  control  gains  are  required  to  be  stored,  because  the  guidance  law  is 
applied  throughout  the  trajectory.  These  correspond  to  gains  associated  with  the 
change  in  the  control  made  at  the  node  where  the  deviation  occurred,  the  gains 
associated  with  the  change  in  the  control  made  at  the  following  node,  and  the 
gains  associated  with  the  change  in  the  time-to-go.  The  result  of  the  gain 
calculations  for  the  change  in  the  control  are  given  in  Table  5.1.  For  a  deviation 
occurring  at  each  node,  the  gains  associated  with  the  change  in  the  control  at  that 
node  arc  shown  in  the  first  row,  while  the  second  row  shows  the  gains  associated 
with  a  change  in  the  control  at  the  following  node.  All  gains  associated  with  Jt 
arc  zero  and  have  been  omitted  from  the  table.  The  gains  associated  with  the 
change  in  the  time-to-go  are  given  in  Table  5.2,  where  again  the  gains  associated 
with  jt  arc  zero  and  have  been  omitted. 

An  obvious  difference  in  the  gains  presented  by  Nowack  [21  and  the  gains 
presented  in  Tables  5. 1-5.2  is  the  number  of  gains.  When  horizontal  velocity  is 
the  reference  variable,  it  is  not  a  state  and  no  deviation  in  horizontal  velocity 
exists  from  the  nominal  trajectory  at  the  update  point.  Therefore,  there  is  no  gain 
associated  with  a  deviation  in  horizontal  velocity.  In  addition,  time  has  been 
included  in  the  state  vector  instead  of  horizontal  velocity.  Because  time  has  no 
impact  on  satisfying  the  constraints,  the  gains  associated  with  time  are  zero  and 
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NODE 

NODE  GAINS 

_ i _ a _ _ * _ 

1 

-0.369E-5 

d.mE-4 

-0.673E-3 

-0.294E-5 

0.539E-4 

-0.575E-3 

2 

-0.184E-5 

0.17iE-4 

-6.454E-$ 

-6.4&E-5 

0.S66E-4 

-0.7$5E-3 

S 

-0.372E-5 

0.265E-4 

-0.510E-3 

-0.617E-5 

0.978E-4 

-0.91  IE-3 

4 

-0.54^-5 

i  0.321E-4 

-0.623E-3 

-0.842E-5 

O.107E-3 

-0.106E-2 

5 

-088&E-5 

0.5216-4 

-0.&06E-3 

-0.126E-4 

0.114E-3 

-0.125E-2 

6 

-0.168E-4 

0.S93E-4 

-0114E-2 

J  -0.174E-4 

0.104E-3  1 

-6.149E-2 

1 

-0.4296-4 

0.178E-3 

-0.194E-2 

-0.1946-4 

0.709E-5 

-0.159E-2 

$ 

-0.248E-3 

0.5&3E-3 

-0.$60E-2 

0.247E-3 

-0.109E-2 

0.271E-2 

Table  5.1  Control  gains  associated  with  a  change  in  the  control 


not  presented  by  Nowack.  On  the  other  hand,  three  gains  are  required  when  time 
is  the  reference  variable.  These  gains  are  required  to  account  for  deviations  from 
the  nominal  trajectory  in  each  of  the  constrained  states.  Only  x  in  the  state  vector 
does  not  impact  the  final  constraints,  and  therefore  only  the  gains  associated  with 
x  are  zero. 

S3  Trajectory  Simulation 

Simulations  using  the  guidance  law  are  compared  with  results  obtained  by 
Nowack  [2].  Tables  5.3  and  5.4  show  results  when  perturbations  in  thrust 
acceleration  or  lunar  gravity  are  added  into  the  system.  Each  table  lists  the  cause 
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NODE 

NODE  GAINS 

y  u  v 

1 

-0.776E-7 

-6.13^4 

i 

-0.8&7E-'? 

-0.556E-4 

r  -b.i41E-4 

i 

-0.l64fe4 

-0.649E-4  “ 

-0.24^E4 

4 

'0.124E-6 

-0.‘>79E-4 

“  -6.2$3e4 

5 

*0.156Ii-<S 

-6.973E-4 

-0.263E-4 

<J 

-6.2f)8E4 

-0.130E-3 

'  -6.iSiE4 

7 

~0.313Ii4 

-6.1<)5E-3 

-0.315E-4 

$ 

I  -0.632E-6 

-0.3S9E-3 

"  -0.418E-4 

Table  5.2  Control  gains  associated  with  time-to-go 


of  the  perturbation  as  well  as  the  optimal  final  time  if  the  cause  of  the  perturbation 
had  been  assumed  when  computing  the  suboptimal  trajectory.  For  each 
perturbation,  the  first  row  shows  the  deviation  in  the  final  values  of  each  of  the 
constrained  states.  The  second  row  shows  the  same  deviations  for  simulations 
accomplished  using  the  guidance  law  presented  by  Nowack  [2], 

The  actual  optimal  final  time,  shown  in  Tables  5.3  and  5.4,  is  found  using 
the  same  optimization  code  that  solved  for  the  nominal  control  and  optimal  final 
time.  The  optimization  code  determines  the  actual  final  time  when  the  nominal 
value  for  thrust  acceleration  or  lunar  gravity  is  replaced  by  the  perturbed  value  for 
thrust  acceleration  or  lunar  gravity.  The  perturbed  final  times  of  both  guidance 
laws  are  close  to  the  actual  optimal  final  time.  The  difference  in  the  perturbed 
final  time  and  the  actual  final  time  is  under  0. 1  seconds  in  all  cases. 

The  results  of  the  simulation  also  show  no  large  disparity  between 
guidance  laws  in  terms  of  meeting  endpoint  constraints.  Both  laws  handle 
deviations  caused  by  thrust  acceleration  perturbations  better,  as  evidenced  by 
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smaller  position  errors.  Velocity  errors  are  approximately  the  same  regardless  of 
the  cause  of  the  deviation  or  the  guidance  law  used. 

A  comparison  of  the  guidance  laws  throughout  the  trajectory  is  presented 
in  Figs.  5. 1-5.8.  Figures  5.1  and  5.4  show  the  deviation  in  vertical  velocity  of  the 
perturbed  trajectory  from  the  nominal  trajectory  when  the  control  updates  are 
made  for  each  of  the  guidance  laws.  The  perturbed  trajectory  remains  closer  to 
the  nominal  trajectory  with  time  as  the  reference  variable,  regardless  of  the 
perturbation  in  the  system.  The  deviations  in  vertical  velocity,  shown  in  Figs.  5.2 
and  5.5,  show  only  a  slight  improvement  with  time  as  the  reference  variable  and 
when  the  deviations  are  caused  by  perturbations  in  lunar  gravity.  The  maximum 
deviation  in  vertical  velocity  caused  by  errors  in  thrust  acceleration  are  slightly 
reduced  with  time  as  the  reference  variable.  However,  deviations  from  the 
nominal  trajectory  are  not  always  smaller  with  time  as  the  reference  variable  than 
when  horizontal  velocity  is  the  reference  variable. 

The  deviation  in  horizontal  velocity  from  the  nominal  trajectory  is  zero 
when  horizontal  velocity  is  the  reference  variable.  Therefore,  it  is  not  shown  in 
Figs.  5.3  and  5.7.  The  magnitude  of  the  deviations  in  horizontal  velocity  with 
time  as  the  reference  variable,  however,  is  smaller  than  deviations  in  vertical 
velocity. 

Figures  5.4  and  5.8  show  the  norm  of  the  error  of  the  three  constrained 
states  throughout  the  trajectory.  Both  guidance  laws  have  similar  norms 
throughout  the  trajectory,  regardless  of  the  cause  of  the  deviation.  The  norm  with 
time  as  the  reference  variable  is  slightly  smaller  throughout  the  trajectory  when 
the  deviations  are  caused  by  perturbations  in  lunar  gravity.  Perturbations  in  thrust 
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accelerations  lead  to  smaller  norms  overall.  However,  neither  guidance  law  yields 
a  smaller  norm  throughout  the  trajectory.  Both  guidance  laws  show  an  obvious 
increase  in  the  norm  when  the  last  control  node  is  passed. 


a 


Reference 


Perturbed 


Final  Deviation  in  States 


fc?) 


19.760 


19.968 


20.176 


20.384 


20.592 


21.008 


21.216 


21.424 


21.632 


21.840 


Variable  Final  Time 
(sec) 


287.957 

8.958 

0.241 

-2.868 

287.947 

10.530 

0.0 

-2.821 

284.764 

7.724 

0.192 

-2.296 

284.756 

8.931 

0.0 

-2.263 

281.645 

6.209 

0.143 

-1.716 

281.639 

7.086 

0.0 

-1.703 

278.596 

4.445 

0.096 

-1.152 

278.593 

4.988 

0.0 

-1.140 

275.619 

2.390 

0.049 

-0.588 

275.616 

2.630 

0.0 

-0.573 

269.857 

-2.664 

-0.052 

0.588 

269.860 

-2.909 

0.0 

0.579 

267.070 

-5.695 

-0.124 

1.200 

267.076 

-6.111 

0.0 

1.167 

264.343 

-9.078 

-0.155 

1.824 

264.351 

-9.615 

0.0 

2.369 

261.672 

-12.783 

-0.254 

2.434 

261.684 

-13.435 

0.0 

2.369 

259.053 

-16.990 

-0.418 

3.127 

259.073 

-17.582 

0.0 

2.987 

Table  5.3  Results  for  errors  in  modeling  thrust  acceleration 


5.054 


5.107 


5.160 


5.214 


5.267 


5.373 


5.426 


5.480 


5.533 


5.586 
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Table  5.4  Results  for  errors  in  modeling  lunar  gravity 
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Figure  5.1  Deviation  in  vertical  position  for  a  thrust  acceleration  perturi 


Figure  5.2  Deviation  in  vertical  velocity  for  a  thrust  acceleration  perturbation 
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Figure  5.7  Deviation  in  horizontal  velocity  for  a  lunar  gravity  perturbation 


Figure  5.8  Norm  of  error  for  a  lunar  gravity  perturbation 
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Chapter  6 


Conclusions 

The  purpose  of  this  research  is  to  demonstrate  the  neighboring  extremal 
guidance  law  using  time  as  the  reference  variable.  The  use  of  horizontal  velocity 
as  the  reference  variable  has  been  previously  demonstrated  by  Nowack  [2].  By 
using  current  time  as  the  reference  variable,  the  location  of  the  vehicle  on  the 
perturbed  trajectory  is  unknown,  since  the  final  time  is  continually  changing. 
Therefore,  a  one-dimensional  search  is  used  to  locate  the  vehicle  on  the  perturbed 
trajectory  by  relating  time-to-go  on  the  perturbed  trajectory  to  the  same  time-to-go 
on  the  nominal  trajectory.  This  provides  the  appropriate  reference  time  to 
determine  the  nominal  control  and  the  change  in  the  control. 

The  guidance  law,  with  time  as  the  reference  variable,  is  demonstrated 
using  the  lunar  launch  problem.  The  lunar  launch  problem  involves  placing  a 
vehicle  into  lunar  orbit  in  minimum  time  while  insuring  that  certain  orbit  entry 
parameters  are  satisfied.  The  perturbed  trajectory  is  integrated  using  equations  of 
motion  that  have  been  altered  from  those  equations  of  motion  used  to  generate  the 
nominal  trajectory.  These  altered  equations  of  motion  are  determined  by 
changing  the  values  of  thrust  acceleration  or  lunar  gravity.  The  states  at  each 
sample  time  are  then  compared  to  the  states  on  the  nominal  trajectory.  A 
reference  time  is  found  on  the  nominal  trajectory  such  that  the  time-to-go  on  both 
trajectories  is  the  same.  The  nominal  control  at  this  reference  time  and  the  change 
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in  the  nominal  control  caused  by  deviations  in  the  states  are  used  to  determine  the 
new  control  history. 

The  results  of  this  simulation  are  compared  with  results  obtained  using 
horizontal  velocity  as  the  reference  variable.  In  all  cases,  the  results  are 
comparable.  When  thrust  acceleration  is  perturbed  and  time  is  used  as  the 
reference  variable,  the  final  value  of  the  states  differ  from  Nowack's  [2]  results  by 
less  then  2  feet  in  position  and  by  less  the  0.5  ft/sec  in  velocity.  The  difierences 
are  slightly  larger,  10  feet  and  1  ft/sec,  when  the  perturbations  are  caused  by 
errors  in  lunar  gravity. 

Both  laws  display  a  noticeable  growth  in  the  norm  of  the  error  after  the 
last  node  point  is  reached.  This  is  expected,  because  the  change  in  the  control 
over  the  last  interval  is  an  extrapolation  from  the  previous  interval.  This  leads  to 
an  obvious  conclusion  that  meeting  the  final  constraints  exactly  is  an 
impossibility.  With  this  is  mind,  no  definition  of  an  acceptable  error  has  been 
given,  making  it  difficult  to  determine  how  well  the  guidance  law  behaves.  One 
definition,  used  by  Helffich  [3],  has  the  criterion  that  the  percentage  error  in  the 
final  states  should  not  exceed  the  percentage  error  in  the  perturbations.  This 
criteria  seems  too  restrictive  when  the  final  constraint  is  zero.  Therefore,  realistic 
problems  need  to  be  used  in  which  realistic  and  acceptable  errors  are  defined. 

Areas  of  further  research  include  understanding  time-to-go  and  its 
relationship  to  interpreting  the  change  in  the  control.  This  research  has  tacitly 
assumed,  and  not  satisfactorily  explained,  a  relationship  between  time-to-go  and 
the  determination  of  the  appropriate  change  in  the  control.  However,  it  has 
demonstrated  through  results  that  a  relationship  exists.  Also,  an  investigation  of 


better  terminal  guidance  is  required.  Either  increasing  the  number  of  control 
nodes  or  shifting  the  nodes  closer  to  the  final  constraints  may  increase  the 
accuracy  of  the  terminal  phase  of  the  trajectory.  This  research  has  used  linear 
extrapolation  over  the  last  interval,  but  it  is  known  that  the  gains  increase  to 
infinity  as  the  final  constraints  are  reached.  Therefore,  a  different  method  of 
extrapolating  the  gains  over  the  last  interval  might  yield  better  results. 
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